BeeResp_15Feb
Callin Switzer
February 15, 2017; Update 1 March 2017; Update March 16
Summary:
Two avenues of analysis:
- First we see how \(\Delta\)load affects \(\Delta\)ArcLength^2 and \(\Delta\)frequency^2, accounting for bee size and treatment order
- We check to see how \(\Delta\)frequency and \(\Delta\)ArcLength are associated with \(\Delta\)MetabolicRate, while accounting for \(\Delta\)load
Measured variables:
- Bee ID: Name given to an invididual (subscript in equations)
- Order: This is the trial number per individual, either 1 or 2
- Treatment: Either light or heavy- H is with both nectar and external load, L is just nectar load (subscript in equations)
- MStarved: The empty mass of the bee after being starved for days, until it no longer buzzes when prodded in grams
- M2: Total mass of the bee and load at beginning of trial in grams
- MF: The mass of the bee and load at end of trial in grams
- ITspan=intertegular span in m
- S=area of both forewings in \(m^2\)
- MetR= Metabolic rate in mL \(CO_2\) \(hr^{-1}\)
- Freq= wing beat frequency in Hz
- Amp= stroke amplitude of forewing in degrees
- L= length of forewing from wing base to tip in m
Calculated variables:
- MT: Mtrial, total mass, =(M2+MF)/2 = the average mass of the course of the trial
- Load= Mtrial- Mstarved
- Perload= Percent load, %load (Load/Mstarved)x100
- ArcL= (.75 \(\cdotp\) wlength) (Amp (\(\pi\)/180)) pi/180 converts to radians
- U= wing velocity, arcL\(\cdotp\) frq \(\cdotp\) 2 (the 2 comes from the fact that frequency consists of both an up and down stroke, vel=∆x/∆t, ∆x=arclength, ∆t=1/(freg \(\cdotp\) 2)
- Freqsq=\(Freq^2\)
- ArcLsq=\(ArcL^2\)
- Force=\(U^2 \cdotp S\)
Install packages and read data
ipak <- function(pkg){
new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
sapply(pkg, require, character.only = TRUE)
}
packages <- c('ggplot2', 'car', 'lme4', 'gsheet', "MASS",'reshape2',
'influence.ME', 'sjPlot', "effects", 'visreg', 'viridis')
ipak(packages)## ggplot2 car lme4 gsheet MASS
## TRUE TRUE TRUE TRUE TRUE
## reshape2 influence.ME sjPlot effects visreg
## TRUE TRUE TRUE TRUE TRUE
## viridis
## TRUE# read in data -- google sheet called "Final Bee Resp Data"
url = 'https://docs.google.com/spreadsheets/d/1wT-QxSJJElhhJcIXlg2hDpFKHbLyFGuqNNj2iYvA8Vo/edit?usp=sharing'
bdta <- data.frame(gsheet2tbl(url))
summary(bdta)## BeeID order Treatment Mstarved
## Length:60 Min. :1.0 Length:60 Min. :0.0767
## Class :character 1st Qu.:1.0 Class :character 1st Qu.:0.1213
## Mode :character Median :1.5 Mode :character Median :0.1378
## Mean :1.5 Mean :0.1411
## 3rd Qu.:2.0 3rd Qu.:0.1670
## Max. :2.0 Max. :0.1909
## M2 MF Itspan S
## Min. :0.1062 Min. :0.0969 Min. :0.003860 Min. :4.010e-05
## 1st Qu.:0.1808 1st Qu.:0.1762 1st Qu.:0.004400 1st Qu.:5.590e-05
## Median :0.2191 Median :0.2107 Median :0.004620 Median :6.130e-05
## Mean :0.2245 Mean :0.2154 Mean :0.004594 Mean :6.141e-05
## 3rd Qu.:0.2646 3rd Qu.:0.2576 3rd Qu.:0.004870 3rd Qu.:6.750e-05
## Max. :0.3550 Max. :0.3440 Max. :0.005180 Max. :7.810e-05
## MetR freq amp L
## Min. : 5.507 Min. :163.9 Min. : 96.88 Min. :0.008691
## 1st Qu.: 8.935 1st Qu.:173.0 1st Qu.:114.16 1st Qu.:0.010570
## Median :10.666 Median :179.1 Median :125.96 Median :0.011017
## Mean :10.600 Mean :178.9 Mean :123.31 Mean :0.010911
## 3rd Qu.:12.297 3rd Qu.:185.3 3rd Qu.:133.61 3rd Qu.:0.011444
## Max. :16.475 Max. :194.7 Max. :144.78 Max. :0.012118Calculate variables
bdta <- within(bdta, {
MT = (M2 + MF)/2
load = MT - Mstarved
perLoad = load / Mstarved * 100
arcL = (0.75 * L) * (amp * (pi / 180))
U = arcL * freq * 2
freq2 = freq^2
arcL2 = arcL^2
frce = U^2 * S
# convert to factor variables
order = as.factor(as.character(order))
Treatment = as.factor(as.character(Treatment))
BeeID = as.factor(as.character(BeeID))
})Create a new dataframe that calculates the changes for each bee
# convert from long to wide
newDF <- data.frame()
colsTocalc = c("order", "M2", "MF", "MetR", "freq", "amp", "load", "MT", "perLoad", "frce", "arcL2", "freq2", "U", "arcL")
for(varb in colsTocalc){
data_wide <- dcast(bdta, BeeID + Mstarved + S + Itspan + L ~
Treatment, value.var=c(varb))
colnames(data_wide)[6:7] = paste(varb, colnames(data_wide)[6:7], sep = "_")
if(varb == colsTocalc[1]){
newDF <- data_wide
}
else newDF <- merge(newDF, data_wide, all.y = TRUE)
}
head(newDF)## BeeID Mstarved S Itspan L order_H order_L M2_H M2_L
## 1 E01 0.1577 6.46e-05 0.00464 0.01116843 2 1 0.2634 0.1895
## 2 E03 0.1732 6.92e-05 0.00487 0.01144364 1 2 0.3106 0.2105
## 3 E05 0.1755 7.06e-05 0.00497 0.01168730 2 1 0.2776 0.2006
## 4 E09 0.1135 4.88e-05 0.00440 0.00986884 1 2 0.2282 0.1425
## 5 E10 0.1269 5.91e-05 0.00462 0.01060036 1 2 0.2166 0.1460
## 6 E11 0.1350 5.88e-05 0.00464 0.01070018 1 2 0.2223 0.1599
## MF_H MF_L MetR_H MetR_L freq_H freq_L amp_H amp_L
## 1 0.2559 0.1812 12.54241 10.393483 186.9973 179.3308 128.7145 111.49875
## 2 0.3079 0.2023 11.78387 9.214622 174.9374 167.8412 144.7784 119.98953
## 3 0.2742 0.1975 12.20102 8.981412 178.8896 166.2385 125.6778 116.22053
## 4 0.2084 0.1390 9.47114 7.581104 191.8348 190.3535 139.3246 103.30703
## 5 0.2099 0.1422 11.03371 7.159840 186.1762 169.5451 128.8257 96.88003
## 6 0.2179 0.1579 10.23580 6.806696 183.0259 168.2223 128.7916 116.18262
## load_H load_L MT_H MT_L perLoad_H perLoad_L frce_H
## 1 0.10195 0.02765 0.25965 0.18535 64.64807 17.53329 0.003199482
## 2 0.13605 0.03320 0.30925 0.20640 78.55081 19.16859 0.003984230
## 3 0.10040 0.02355 0.27590 0.19905 57.20798 13.41880 0.003340857
## 4 0.10480 0.02725 0.21830 0.14075 92.33480 24.00881 0.002327019
## 5 0.08635 0.01720 0.21325 0.14410 68.04571 13.55398 0.002618299
## 6 0.08510 0.02390 0.22010 0.15890 63.03704 17.70370 0.002563877
## frce_L arcL2_H arcL2_L freq2_H freq2_L U_H
## 1 0.002208025 0.0003540922 0.0002657061 34967.99 32159.55 7.037583
## 2 0.002519159 0.0004703414 0.0003230668 30603.08 28170.67 7.587858
## 3 0.002467172 0.0003696773 0.0003161342 32001.50 27635.24 6.879020
## 4 0.001259710 0.0003239404 0.0001781022 36800.61 36234.44 6.905419
## 5 0.001228018 0.0003195386 0.0001807120 34661.58 28745.53 6.656039
## 6 0.001762569 0.0003254128 0.0002648146 33498.49 28298.73 6.603283
## U_L arcL_H arcL_L
## 1 5.846363 0.01881734 0.01630050
## 2 6.033576 0.02168736 0.01797406
## 3 5.911496 0.01922699 0.01778016
## 4 5.080722 0.01799834 0.01334549
## 5 4.558361 0.01787564 0.01344292
## 6 5.475004 0.01803920 0.01627313Calculate \(\Delta\) variables
newDF <- within(newDF, {
deltaPercLoad = perLoad_H - perLoad_L
avgPercLoad = (perLoad_H + perLoad_L) / 2
deltaMetR = MetR_H - MetR_L
deltaFrq2 = freq2_H - freq2_L
deltaArcL = arcL_H - arcL_L
deltaArcL2 = arcL2_H - arcL2_L
deltaFreq2Perc = deltaFrq2 / deltaPercLoad
deltaLoad = scale(load_H - load_L, center = TRUE, scale = FALSE)
dLoad_nonCent = load_H - load_L
deltaLoad2 <- deltaLoad^2
})
plot(newDF$deltaArcL2 ~ newDF$deltaFrq2)
# compare with Susie's calculations
b2 <- read.csv("~/Desktop/b2.csv")
b2 <- b2[!is.na(b2$deltaperload), ]
b2 <- b2[order(b2$BeeID, b2$Treatment), ]
plot(b2$AverageperLoad, newDF$avgPercLoad)
abline(0,1)
plot(b2$deltaload, newDF$deltaLoad)
abline(0,1)
Use PCA to combine bee sizes into a single predictor
# principle components
aa = prcomp(newDF[, c("Mstarved", "Itspan", "S", "L")], center = TRUE, scale = TRUE)
summary(aa) # 1st pc explains ~95% of the variance in the three measurements of size## Importance of components:
## PC1 PC2 PC3 PC4
## Standard deviation 1.9422 0.37234 0.22746 0.19332
## Proportion of Variance 0.9431 0.03466 0.01293 0.00934
## Cumulative Proportion 0.9431 0.97772 0.99066 1.00000biplot(aa) # shows that all three size measurement are correlated
# note, I changed the signs of the predictions so that higher PC1 values
# correspond to bigger bees
p1 = -predict(aa)[,1]
# add PC1 scores to dataset
newDF$size_pc1 = p1
newDF$size_pc1_2 = newDF$size_pc1^2
# show scatterplot matrix to see correlations among size predictors
car::scatterplotMatrix(newDF[, c("Mstarved", "Itspan", "S", "L", "size_pc1")])
1. First we see how \(\Delta\)load affects \(\Delta\)ArcLength accounting for bee size and treatment order
This is the model selection procedure that was used: 1. Fit a large model with all two-way interactions and squared terms 2. Remove non-significant predictors, starting with interactions, squared terms, and then main effects
# reformat order so that it is more interpretable
library(plyr)
newDF$order_1 <- mapvalues(newDF$order_H, from = c(2, 1), to = c("loadedSecond", "loadedFirst"))
# fit full model
m1 <- lm(deltaArcL2 ~ (deltaLoad + size_pc1 + order_1)^2 +
size_pc1_2 + deltaLoad2, data = newDF)
summary(m1)##
## Call:
## lm(formula = deltaArcL2 ~ (deltaLoad + size_pc1 + order_1)^2 +
## size_pc1_2 + deltaLoad2, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.326e-05 -1.951e-05 3.244e-06 1.635e-05 4.401e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.818e-05 1.300e-05 7.551 2.05e-07 ***
## deltaLoad 1.859e-03 1.178e-03 1.578 0.130
## size_pc1 -1.804e-05 1.263e-05 -1.428 0.168
## order_1loadedSecond 7.310e-07 1.978e-05 0.037 0.971
## size_pc1_2 -3.660e-06 4.905e-06 -0.746 0.464
## deltaLoad2 3.807e-02 4.366e-02 0.872 0.393
## deltaLoad:size_pc1 -2.300e-06 8.387e-04 -0.003 0.998
## deltaLoad:order_1loadedSecond -6.487e-04 2.225e-03 -0.292 0.774
## size_pc1:order_1loadedSecond 1.204e-05 2.099e-05 0.574 0.572
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.03e-05 on 21 degrees of freedom
## Multiple R-squared: 0.575, Adjusted R-squared: 0.413
## F-statistic: 3.551 on 8 and 21 DF, p-value: 0.009387m2 <- update(m1, .~. - deltaLoad:size_pc1)
anova(m1, m2)## Analysis of Variance Table
##
## Model 1: deltaArcL2 ~ (deltaLoad + size_pc1 + order_1)^2 + size_pc1_2 +
## deltaLoad2
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 +
## deltaLoad:order_1 + size_pc1:order_1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 21 1.9281e-08
## 2 22 1.9281e-08 -1 -6.9071e-15 0 0.9978summary(m2)##
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 +
## deltaLoad2 + deltaLoad:order_1 + size_pc1:order_1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.327e-05 -1.953e-05 3.263e-06 1.634e-05 4.402e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.817e-05 1.252e-05 7.840 8.24e-08 ***
## deltaLoad 1.860e-03 9.935e-04 1.872 0.0745 .
## size_pc1 -1.806e-05 7.886e-06 -2.290 0.0320 *
## order_1loadedSecond 7.471e-07 1.846e-05 0.040 0.9681
## size_pc1_2 -3.672e-06 1.745e-06 -2.104 0.0470 *
## deltaLoad2 3.800e-02 3.336e-02 1.139 0.2669
## deltaLoad:order_1loadedSecond -6.509e-04 2.022e-03 -0.322 0.7505
## size_pc1:order_1loadedSecond 1.209e-05 1.187e-05 1.019 0.3194
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.96e-05 on 22 degrees of freedom
## Multiple R-squared: 0.575, Adjusted R-squared: 0.4397
## F-statistic: 4.251 on 7 and 22 DF, p-value: 0.004168m3 <- update(m2, .~. - deltaLoad:order_1)
anova(m2, m3)## Analysis of Variance Table
##
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 +
## deltaLoad:order_1 + size_pc1:order_1
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 +
## size_pc1:order_1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 22 1.9281e-08
## 2 23 1.9371e-08 -1 -9.0841e-11 0.1037 0.7505summary(m3)##
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 +
## deltaLoad2 + size_pc1:order_1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.414e-05 -1.856e-05 3.421e-06 1.715e-05 4.272e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.768e-05 1.218e-05 8.017 4.13e-08 ***
## deltaLoad 1.625e-03 6.585e-04 2.467 0.0215 *
## size_pc1 -1.728e-05 7.355e-06 -2.349 0.0278 *
## order_1loadedSecond 3.683e-06 1.573e-05 0.234 0.8170
## size_pc1_2 -3.816e-06 1.654e-06 -2.308 0.0304 *
## deltaLoad2 4.504e-02 2.467e-02 1.825 0.0809 .
## size_pc1:order_1loadedSecond 1.003e-05 9.807e-06 1.023 0.3169
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.902e-05 on 23 degrees of freedom
## Multiple R-squared: 0.573, Adjusted R-squared: 0.4615
## F-statistic: 5.143 on 6 and 23 DF, p-value: 0.001756m4 <- update(m3, .~. - size_pc1:order_1)
anova(m3, m4) ## Analysis of Variance Table
##
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 +
## size_pc1:order_1
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 23 1.9371e-08
## 2 24 2.0253e-08 -1 -8.8175e-10 1.0469 0.3169summary(m4)##
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 +
## deltaLoad2, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.427e-05 -1.874e-05 -5.000e-09 1.783e-05 4.254e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.024e-04 1.128e-05 9.079 3.13e-09 ***
## deltaLoad 1.639e-03 6.590e-04 2.486 0.0203 *
## size_pc1 -1.180e-05 5.042e-06 -2.340 0.0279 *
## order_1loadedSecond -5.411e-07 1.519e-05 -0.036 0.9719
## size_pc1_2 -2.627e-06 1.178e-06 -2.231 0.0353 *
## deltaLoad2 2.613e-02 1.636e-02 1.597 0.1233
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.905e-05 on 24 degrees of freedom
## Multiple R-squared: 0.5535, Adjusted R-squared: 0.4605
## F-statistic: 5.951 on 5 and 24 DF, p-value: 0.001027m5 <- update(m4, .~. - deltaLoad2)
anova(m4,m5) ## Analysis of Variance Table
##
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 24 2.0253e-08
## 2 25 2.2406e-08 -1 -2.1529e-09 2.5511 0.1233summary(m5)##
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2,
## data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.782e-05 -1.748e-05 9.200e-07 1.678e-05 4.278e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.074e-04 1.117e-05 9.623 6.93e-10 ***
## deltaLoad 2.095e-03 6.122e-04 3.422 0.00215 **
## size_pc1 -1.388e-05 5.018e-06 -2.767 0.01050 *
## order_1loadedSecond 2.783e-07 1.565e-05 0.018 0.98595
## size_pc1_2 -2.037e-06 1.152e-06 -1.768 0.08930 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.994e-05 on 25 degrees of freedom
## Multiple R-squared: 0.5061, Adjusted R-squared: 0.427
## F-statistic: 6.403 on 4 and 25 DF, p-value: 0.001086m6 <- update(m5, .~. - size_pc1_2)
anova(m6, m5) ## Analysis of Variance Table
##
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1 + order_1
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 26 2.5207e-08
## 2 25 2.2406e-08 1 2.8009e-09 3.1251 0.0893 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m6)##
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1 + order_1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.154e-05 -1.781e-05 -4.100e-07 1.849e-05 4.983e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.844e-05 1.034e-05 9.525 5.79e-10 ***
## deltaLoad 2.148e-03 6.359e-04 3.378 0.00231 **
## size_pc1 -1.321e-05 5.205e-06 -2.539 0.01745 *
## order_1loadedSecond 3.224e-06 1.618e-05 0.199 0.84364
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.114e-05 on 26 degrees of freedom
## Multiple R-squared: 0.4443, Adjusted R-squared: 0.3802
## F-statistic: 6.929 on 3 and 26 DF, p-value: 0.001402m7 <- update(m6, .~. - order_1)
anova(m7, m6)## Analysis of Variance Table
##
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 27 2.5245e-08
## 2 26 2.5207e-08 1 3.8477e-11 0.0397 0.8436summary(m7)##
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.041e-05 -1.877e-05 -2.970e-07 1.831e-05 4.778e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.002e-04 5.583e-06 17.940 < 2e-16 ***
## deltaLoad 2.059e-03 4.439e-04 4.638 8.06e-05 ***
## size_pc1 -1.256e-05 3.967e-06 -3.166 0.00381 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.058e-05 on 27 degrees of freedom
## Multiple R-squared: 0.4435, Adjusted R-squared: 0.4022
## F-statistic: 10.76 on 2 and 27 DF, p-value: 0.0003666anova(m7, update(m7, .~. - deltaLoad)) # p-value for delta load## Analysis of Variance Table
##
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1
## Model 2: deltaArcL2 ~ size_pc1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 27 2.5245e-08
## 2 28 4.5359e-08 -1 -2.0114e-08 21.512 8.055e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m7, update(m7, .~. - size_pc1)) # p-value for size## Analysis of Variance Table
##
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1
## Model 2: deltaArcL2 ~ deltaLoad
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 27 2.5245e-08
## 2 28 3.4619e-08 -1 -9.3737e-09 10.025 0.003808 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m7)##
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.041e-05 -1.877e-05 -2.970e-07 1.831e-05 4.778e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.002e-04 5.583e-06 17.940 < 2e-16 ***
## deltaLoad 2.059e-03 4.439e-04 4.638 8.06e-05 ***
## size_pc1 -1.256e-05 3.967e-06 -3.166 0.00381 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.058e-05 on 27 degrees of freedom
## Multiple R-squared: 0.4435, Adjusted R-squared: 0.4022
## F-statistic: 10.76 on 2 and 27 DF, p-value: 0.0003666# refit model with non-centered version of deltaLoad
# this model will have a different intercept, but same p-values for slopes
summary(update(m7, .~. - deltaLoad + dLoad_nonCent)) # final model for paper##
## Call:
## lm(formula = deltaArcL2 ~ size_pc1 + dLoad_nonCent, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.041e-05 -1.877e-05 -2.970e-07 1.831e-05 4.778e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.822e-05 3.460e-05 -1.683 0.10397
## size_pc1 -1.256e-05 3.967e-06 -3.166 0.00381 **
## dLoad_nonCent 2.059e-03 4.439e-04 4.638 8.06e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.058e-05 on 27 degrees of freedom
## Multiple R-squared: 0.4435, Adjusted R-squared: 0.4022
## F-statistic: 10.76 on 2 and 27 DF, p-value: 0.0003666model diagnostics
par(mfrow = c(2,2))
plot(m7, which = 1:4) # no glaring violations, though there are a few fairly influential observations
par(mfrow = c(1,1))
car::vif(m7) ## deltaLoad size_pc1
## 1.841066 1.841066Model visualization
summary(m7)##
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.041e-05 -1.877e-05 -2.970e-07 1.831e-05 4.778e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.002e-04 5.583e-06 17.940 < 2e-16 ***
## deltaLoad 2.059e-03 4.439e-04 4.638 8.06e-05 ***
## size_pc1 -1.256e-05 3.967e-06 -3.166 0.00381 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.058e-05 on 27 degrees of freedom
## Multiple R-squared: 0.4435, Adjusted R-squared: 0.4022
## F-statistic: 10.76 on 2 and 27 DF, p-value: 0.0003666# calculate partial residuals for deltaLoad
# these are the residuals, minus the effect of detlaLoad
y <- residuals(m7, type = 'partial')[, "deltaLoad"]
# plot partial residuals with base R plotting
plot(x = newDF$deltaLoad, y = y)
# double check to make sure the slope for partial residual plots are the
# same as in the original regression
summary(lm(y ~ newDF$dLoad_nonCent)) ##
## Call:
## lm(formula = y ~ newDF$dLoad_nonCent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.041e-05 -1.877e-05 -2.970e-07 1.831e-05 4.778e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.584e-04 2.531e-05 -6.257 9.20e-07 ***
## newDF$dLoad_nonCent 2.059e-03 3.213e-04 6.409 6.14e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.003e-05 on 28 degrees of freedom
## Multiple R-squared: 0.5946, Adjusted R-squared: 0.5801
## F-statistic: 41.07 on 1 and 28 DF, p-value: 6.137e-07# this is what the raw data look like
plot(x = newDF$deltaLoad, y = newDF$deltaFrq2)
# this package plots the partial residuals
crPlot(m7, variable = "deltaLoad")
# plot with ggplot2
theme_set(theme_classic())
# plot raw data w/ ggplot
ggplot(newDF, aes(x= size_pc1, y = deltaArcL, color = deltaLoad)) +
geom_point()
ggplot(newDF, aes(x= deltaLoad, y = deltaArcL, color = size_pc1)) +
geom_point()
# y axis isn't easily interpretable
ggplot(newDF, aes(x= dLoad_nonCent, y = y)) +
geom_point() +
labs(x = "delta load", y = "partial residuals for delta load \n i.e. delta load effect on arcL^2 \n while subtracing affect of bee size") +
stat_smooth(method = 'lm', se = FALSE)
partialResidDeltaLoad <- data.frame(deltaLoad = newDF$dLoad_nonCent, partResArcL2 = y )Takeaways from model 1:
- A higher deltaload (i.e. a relatively larger “high” load) causes bees to have a larger change in arclength squared. Another way to say this is that increasing the change in load by 1 gram causes an increase in the change in arclength^2 by 0.002059 degrees.
- The larger the bee, the lower the change in arcLength^2 (association, rather than causation), while holding the change in load constant.
- We do not have enough evidence to notice any non-linear relationships in this model.
See how \(\Delta\)load affects \(\Delta\)frequency, accounting for bee size and treatment order
# fit full model
m1 <- lm(deltaFrq2 ~ (deltaLoad + size_pc1 + order_1)^2 +
size_pc1_2 + deltaLoad2, data = newDF)
summary(m1)##
## Call:
## lm(formula = deltaFrq2 ~ (deltaLoad + size_pc1 + order_1)^2 +
## size_pc1_2 + deltaLoad2, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3981.4 -1255.9 -170.3 1165.6 3552.9
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4195.97 909.39 4.614 0.00015 ***
## deltaLoad -106368.80 82390.69 -1.291 0.21073
## size_pc1 666.06 883.26 0.754 0.45917
## order_1loadedSecond -4376.19 1383.48 -3.163 0.00469 **
## size_pc1_2 -120.74 343.06 -0.352 0.72839
## deltaLoad2 -469149.49 3053199.27 -0.154 0.87935
## deltaLoad:size_pc1 6902.45 58657.40 0.118 0.90744
## deltaLoad:order_1loadedSecond 62840.01 155625.35 0.404 0.69045
## size_pc1:order_1loadedSecond 54.39 1467.95 0.037 0.97079
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2119 on 21 degrees of freedom
## Multiple R-squared: 0.5048, Adjusted R-squared: 0.3162
## F-statistic: 2.676 on 8 and 21 DF, p-value: 0.03372car::vif(m1) # some serious multicollinearity## deltaLoad size_pc1 order_1
## 13.203325 19.003708 3.182227
## size_pc1_2 deltaLoad2 deltaLoad:size_pc1
## 18.740891 10.117446 25.115087
## deltaLoad:order_1 size_pc1:order_1
## 11.384208 23.301468m2 <- update(m1, .~. - size_pc1:order_1)
anova(m1, m2)## Analysis of Variance Table
##
## Model 1: deltaFrq2 ~ (deltaLoad + size_pc1 + order_1)^2 + size_pc1_2 +
## deltaLoad2
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 +
## deltaLoad:size_pc1 + deltaLoad:order_1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 21 94310333
## 2 22 94316499 -1 -6165.2 0.0014 0.9708summary(m2)##
## Call:
## lm(formula = deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 +
## deltaLoad2 + deltaLoad:size_pc1 + deltaLoad:order_1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3973.1 -1259.6 -176.2 1149.9 3556.7
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4205.4 852.7 4.932 6.21e-05 ***
## deltaLoad -108212.1 64166.1 -1.686 0.10584
## size_pc1 695.7 367.7 1.892 0.07176 .
## order_1loadedSecond -4386.5 1324.0 -3.313 0.00316 **
## size_pc1_2 -109.9 176.0 -0.625 0.53864
## deltaLoad2 -417980.4 2660489.9 -0.157 0.87659
## deltaLoad:size_pc1 5130.1 33170.6 0.155 0.87850
## deltaLoad:order_1loadedSecond 66238.1 122843.3 0.539 0.59516
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2071 on 22 degrees of freedom
## Multiple R-squared: 0.5048, Adjusted R-squared: 0.3472
## F-statistic: 3.203 on 7 and 22 DF, p-value: 0.01701m3 <- update(m2, .~. - deltaLoad:size_pc1)
anova(m2, m3)## Analysis of Variance Table
##
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 +
## deltaLoad:size_pc1 + deltaLoad:order_1
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 +
## deltaLoad:order_1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 22 94316499
## 2 23 94419044 -1 -102545 0.0239 0.8785summary(m3)##
## Call:
## lm(formula = deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 +
## deltaLoad2 + deltaLoad:order_1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3910.3 -1244.0 -168.5 1173.4 3503.8
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4190.99 829.34 5.053 4.09e-05 ***
## deltaLoad -108478.92 62767.13 -1.728 0.09734 .
## size_pc1 702.89 356.93 1.969 0.06108 .
## order_1loadedSecond -4434.40 1259.70 -3.520 0.00184 **
## size_pc1_2 -88.67 107.58 -0.824 0.41829
## deltaLoad2 -215246.26 2265485.45 -0.095 0.92513
## deltaLoad:order_1loadedSecond 61643.10 116639.62 0.528 0.60222
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2026 on 23 degrees of freedom
## Multiple R-squared: 0.5042, Adjusted R-squared: 0.3749
## F-statistic: 3.899 on 6 and 23 DF, p-value: 0.00785m4 <- update(m3, .~. - deltaLoad:order_1)
anova(m3, m4) ## Analysis of Variance Table
##
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 +
## deltaLoad:order_1
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 23 94419044
## 2 24 95565635 -1 -1146591 0.2793 0.6022summary(m4)##
## Call:
## lm(formula = deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 +
## deltaLoad2, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3667.1 -1308.7 -52.9 1317.9 3482.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.330e+03 7.749e+02 5.587 9.48e-06 ***
## deltaLoad -8.589e+04 4.527e+04 -1.897 0.069889 .
## size_pc1 7.352e+02 3.463e+02 2.123 0.044282 *
## order_1loadedSecond -4.794e+03 1.044e+03 -4.594 0.000117 ***
## size_pc1_2 -5.195e+01 8.090e+01 -0.642 0.526837
## deltaLoad2 -1.250e+06 1.124e+06 -1.112 0.277131
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1995 on 24 degrees of freedom
## Multiple R-squared: 0.4982, Adjusted R-squared: 0.3937
## F-statistic: 4.766 on 5 and 24 DF, p-value: 0.00364m5 <- update(m4, .~. - size_pc1_2)
anova(m4,m5) ## Analysis of Variance Table
##
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + deltaLoad2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 24 95565635
## 2 25 97207799 -1 -1642165 0.4124 0.5268summary(m5)##
## Call:
## lm(formula = deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + deltaLoad2,
## data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3591.2 -1388.0 -154.6 1391.9 3336.4
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4166.0 723.2 5.760 5.29e-06 ***
## deltaLoad -80714.2 44020.7 -1.834 0.078654 .
## size_pc1 732.5 342.2 2.141 0.042251 *
## order_1loadedSecond -4719.5 1024.9 -4.605 0.000104 ***
## deltaLoad2 -1475818.2 1054446.2 -1.400 0.173916
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1972 on 25 degrees of freedom
## Multiple R-squared: 0.4896, Adjusted R-squared: 0.4079
## F-statistic: 5.995 on 4 and 25 DF, p-value: 0.00159m6 <- update(m5, .~. - deltaLoad2)
anova(m6, m5) ## Analysis of Variance Table
##
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + deltaLoad2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 26 104824695
## 2 25 97207799 1 7616896 1.9589 0.1739summary(m6)##
## Call:
## lm(formula = deltaFrq2 ~ deltaLoad + size_pc1 + order_1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3191.2 -1541.4 -264.1 1651.1 3486.0
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3735.3 666.5 5.605 6.87e-06 ***
## deltaLoad -105594.7 41007.6 -2.575 0.016065 *
## size_pc1 861.4 335.6 2.566 0.016382 *
## order_1loadedSecond -4717.7 1043.6 -4.521 0.000119 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2008 on 26 degrees of freedom
## Multiple R-squared: 0.4496, Adjusted R-squared: 0.3861
## F-statistic: 7.08 on 3 and 26 DF, p-value: 0.001244m7 <- update(m6, .~. - size_pc1)
anova(m7, m6) # p-values for size## Analysis of Variance Table
##
## Model 1: deltaFrq2 ~ deltaLoad + order_1
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 27 131380648
## 2 26 104824695 1 26555953 6.5868 0.01638 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m8 <- update(m6, .~. - deltaLoad)
anova(m6, m8) # p-value for deltaLoad## Analysis of Variance Table
##
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1
## Model 2: deltaFrq2 ~ size_pc1 + order_1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 26 104824695
## 2 27 131557512 -1 -26732816 6.6306 0.01607 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m9 <- update(m6, .~. - order_1)
anova(m6, m9) # p-value for order## Analysis of Variance Table
##
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 26 104824695
## 2 27 187214161 -1 -82389465 20.435 0.0001192 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# refit final model with non-centered deltaLoad
m10 <- update(m6, .~. - deltaLoad + dLoad_nonCent)
summary(m10) # final model for paper##
## Call:
## lm(formula = deltaFrq2 ~ size_pc1 + order_1 + dLoad_nonCent,
## data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3191.2 -1541.4 -264.1 1651.1 3486.0
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11857.7 3586.6 3.306 0.002766 **
## size_pc1 861.4 335.6 2.566 0.016382 *
## order_1loadedSecond -4717.7 1043.6 -4.521 0.000119 ***
## dLoad_nonCent -105594.7 41007.6 -2.575 0.016065 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2008 on 26 degrees of freedom
## Multiple R-squared: 0.4496, Adjusted R-squared: 0.3861
## F-statistic: 7.08 on 3 and 26 DF, p-value: 0.001244model diagnostics
par(mfrow = c(2,2))
plot(m10, which = 1:4) # no glaring violations
par(mfrow = c(1,1))
car::vif(m10) # vif is a little high## size_pc1 order_1 dLoad_nonCent
## 3.056456 2.017003 3.643394model visualization
# calculate partial residuals for deltaLoad
# these are the residuals, minus the effect of detlaLoad
y <- residuals(m10, type = 'partial')[, "dLoad_nonCent"]
partialResidDeltaLoad <- cbind(partialResidDeltaLoad, y)
colnames(partialResidDeltaLoad)[3] <- "partResFreq2"
# plot partial residuals with base R plotting
plot(x = newDF$dLoad_nonCent, y = y)
# double check to make sure the slope for partial residual plots are the
# same as in the original regression
summary(lm(y ~ newDF$dLoad_nonCent)) ##
## Call:
## lm(formula = y ~ newDF$dLoad_nonCent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3191.2 -1541.4 -264.1 1651.1 3486.0
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8122 1631 4.980 2.93e-05 ***
## newDF$dLoad_nonCent -105595 20702 -5.101 2.11e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1935 on 28 degrees of freedom
## Multiple R-squared: 0.4816, Adjusted R-squared: 0.4631
## F-statistic: 26.02 on 1 and 28 DF, p-value: 2.106e-05# this is what the raw data look like
plot(x = newDF$dLoad_nonCent, y = newDF$deltaFrq2)
# this package plots the partial residuals
crPlot(m10, variable = "dLoad_nonCent")
# plot with ggplot2
# plot raw data w/ ggplot
ggplot(newDF, aes(x= size_pc1, y = deltaFrq2, color = deltaLoad, shape = order_1)) +
geom_point()
ggplot(newDF, aes(x= deltaLoad, y = deltaFrq2, color = size_pc1, shape = order_1)) +
geom_point()
# y axis isn't easily interpretable
ggplot(newDF, aes(x= dLoad_nonCent, y = y)) +
geom_point() +
labs(x = "delta load", y = "partial residuals for delta load \n i.e. delta load effect on freq^2 \n while subtracing affect of bee size and order") +
stat_smooth(method = 'lm', se = FALSE)
# partial residuals for order
y <- residuals(m10, type = 'partial')[, "order_1"]
ggplot(newDF, aes(x= order_1, y = y)) +
geom_boxplot() +
labs(x = "order", y = "partial residuals for order \n i.e. order effect on freq^2 \n while subtracing affect of bee size and load") +
stat_smooth(method = 'lm', se = FALSE)
# holding bee size and delta load constant, if a bee was loaded second, (confusing, huh?) then it would have a much lower freq^2 than if it was loaded first
summary(m10)##
## Call:
## lm(formula = deltaFrq2 ~ size_pc1 + order_1 + dLoad_nonCent,
## data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3191.2 -1541.4 -264.1 1651.1 3486.0
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11857.7 3586.6 3.306 0.002766 **
## size_pc1 861.4 335.6 2.566 0.016382 *
## order_1loadedSecond -4717.7 1043.6 -4.521 0.000119 ***
## dLoad_nonCent -105594.7 41007.6 -2.575 0.016065 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2008 on 26 degrees of freedom
## Multiple R-squared: 0.4496, Adjusted R-squared: 0.3861
## F-statistic: 7.08 on 3 and 26 DF, p-value: 0.001244Takeaways from model 2:
- Order, size, and deltaLoad are associated with a change in frequency^2. We find no evidence of non-linear relationships or interactions.
- Holding size and order constant, we find that a larger deltaLoad causes a decrease in deltaFrequency^2
- Holding other variables constant, an increase in bee size is associated with a larger deltaFrequency^2.
- Holding other variables constant, if the bee was loaded second, then they had a lower deltaFrequency^2 than if they were loaded first.
# make some predictions
predDF <- data.frame(size_pc1 = 0, order_1 = c("loadedSecond", "loadedFirst"), dLoad_nonCent = mean(newDF$dLoad_nonCent))
predDF$pred_deltaF2 <- predict(m10, predDF)
predDF## size_pc1 order_1 dLoad_nonCent pred_deltaF2
## 1 0 loadedSecond 0.07692 -982.314
## 2 0 loadedFirst 0.07692 3735.3372. We check to see how \(\Delta\)frequency and \(\Delta\)ArcLength are associated with \(\Delta\)MetabolicRate, while accounting for \(\Delta\)load
mm1 <- lm(deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad + size_pc1 + deltaLoad2, data = newDF)
car::vif(mm1)## deltaFrq2 deltaArcL2 deltaLoad size_pc1 deltaLoad2
## 1.060059 1.852143 3.844322 2.679591 1.495113summary(mm1)##
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad +
## size_pc1 + deltaLoad2, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.11198 -0.64060 -0.01055 0.42457 2.97948
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.941e+00 6.182e-01 3.140 0.00444 **
## deltaFrq2 3.543e-04 6.888e-05 5.144 2.89e-05 ***
## deltaArcL2 -2.573e+03 5.900e+03 -0.436 0.66659
## deltaLoad 2.368e+01 1.937e+01 1.223 0.23337
## size_pc1 -5.482e-02 1.445e-01 -0.379 0.70776
## deltaLoad2 4.440e+02 5.114e+02 0.868 0.39381
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9233 on 24 degrees of freedom
## Multiple R-squared: 0.5937, Adjusted R-squared: 0.5091
## F-statistic: 7.014 on 5 and 24 DF, p-value: 0.0003629mm2 <- update(mm1, .~. - size_pc1)
anova(mm1, mm2)## Analysis of Variance Table
##
## Model 1: deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad + size_pc1 + deltaLoad2
## Model 2: deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad + deltaLoad2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 24 20.460
## 2 25 20.582 -1 -0.12268 0.1439 0.7078summary(mm2)##
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad +
## deltaLoad2, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.06884 -0.60009 -0.03648 0.36122 3.09536
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.821e+00 5.223e-01 3.487 0.00182 **
## deltaFrq2 3.571e-04 6.730e-05 5.306 1.69e-05 ***
## deltaArcL2 -1.543e+03 5.147e+03 -0.300 0.76682
## deltaLoad 1.788e+01 1.168e+01 1.530 0.13849
## deltaLoad2 4.891e+02 4.888e+02 1.001 0.32665
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9074 on 25 degrees of freedom
## Multiple R-squared: 0.5913, Adjusted R-squared: 0.5259
## F-statistic: 9.041 on 4 and 25 DF, p-value: 0.0001166mm3 <- update(mm2, .~. - deltaArcL2)
anova(mm3, mm2)## Analysis of Variance Table
##
## Model 1: deltaMetR ~ deltaFrq2 + deltaLoad + deltaLoad2
## Model 2: deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad + deltaLoad2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 26 20.656
## 2 25 20.582 1 0.07399 0.0899 0.7668summary(mm3)##
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + deltaLoad + deltaLoad2,
## data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.05850 -0.58135 -0.00991 0.37733 3.10956
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.682e+00 2.327e-01 7.228 1.12e-07 ***
## deltaFrq2 3.560e-04 6.602e-05 5.393 1.19e-05 ***
## deltaLoad 1.667e+01 1.077e+01 1.548 0.134
## deltaLoad2 4.422e+02 4.550e+02 0.972 0.340
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8913 on 26 degrees of freedom
## Multiple R-squared: 0.5898, Adjusted R-squared: 0.5425
## F-statistic: 12.46 on 3 and 26 DF, p-value: 3.065e-05mm4 <- update(mm3, .~. - deltaLoad2)
anova(mm3, mm4)## Analysis of Variance Table
##
## Model 1: deltaMetR ~ deltaFrq2 + deltaLoad + deltaLoad2
## Model 2: deltaMetR ~ deltaFrq2 + deltaLoad
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 26 20.656
## 2 27 21.407 -1 -0.75053 0.9447 0.34summary(mm4)##
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + deltaLoad, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.17958 -0.59123 -0.06454 0.34379 3.00903
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.824e+00 1.808e-01 10.086 1.18e-10 ***
## deltaFrq2 3.451e-04 6.498e-05 5.310 1.32e-05 ***
## deltaLoad 2.140e+01 9.595e+00 2.231 0.0342 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8904 on 27 degrees of freedom
## Multiple R-squared: 0.5749, Adjusted R-squared: 0.5434
## F-statistic: 18.26 on 2 and 27 DF, p-value: 9.654e-06mm5 <- update(mm4, .~. - deltaLoad)
anova(mm5, mm4) # p-value for load## Analysis of Variance Table
##
## Model 1: deltaMetR ~ deltaFrq2
## Model 2: deltaMetR ~ deltaFrq2 + deltaLoad
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 28 25.352
## 2 27 21.407 1 3.9448 4.9754 0.03421 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mm6 <- update(mm4, .~. - deltaFrq2)
anova(mm4, mm6) # p-value for deltafrq2## Analysis of Variance Table
##
## Model 1: deltaMetR ~ deltaFrq2 + deltaLoad
## Model 2: deltaMetR ~ deltaLoad
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 27 21.407
## 2 28 43.766 -1 -22.359 28.201 1.324e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# refit model with non-centered load
mm7 <- update(mm4, .~. - deltaLoad + dLoad_nonCent)
summary(mm7) # final model for paper##
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + dLoad_nonCent, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.17958 -0.59123 -0.06454 0.34379 3.00903
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.777e-01 7.507e-01 0.237 0.8146
## deltaFrq2 3.451e-04 6.498e-05 5.310 1.32e-05 ***
## dLoad_nonCent 2.140e+01 9.595e+00 2.231 0.0342 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8904 on 27 degrees of freedom
## Multiple R-squared: 0.5749, Adjusted R-squared: 0.5434
## F-statistic: 18.26 on 2 and 27 DF, p-value: 9.654e-06# make a partial residual plot for delta arcl2 (not sure if this is the best way)
mm8 <- update(mm7, .~. + deltaArcL2)
summary(mm8)##
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + dLoad_nonCent + deltaArcL2,
## data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.17804 -0.59319 -0.06199 0.34225 3.01069
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.762e-01 7.684e-01 0.229 0.8204
## deltaFrq2 3.451e-04 6.622e-05 5.211 1.93e-05 ***
## dLoad_nonCent 2.129e+01 1.118e+01 1.905 0.0679 .
## deltaArcL2 1.032e+02 4.877e+03 0.021 0.9833
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9074 on 26 degrees of freedom
## Multiple R-squared: 0.5749, Adjusted R-squared: 0.5259
## F-statistic: 11.72 on 3 and 26 DF, p-value: 4.818e-05crPlot(mm8, variable = "deltaArcL2")
crPlot(mm8, variable = "dLoad_nonCent")
crPlot(mm8, variable = "deltaFrq2")
# calculate partial residuals from model 7 for deltaload and deltafreq2
prDload <- residuals(mm7, type = 'partial')[, "dLoad_nonCent"]
prDfreq2 <- residuals(mm7, type = 'partial')[, "deltaFrq2"]
# plot partial residuals with base R plotting
plot(x = newDF$dLoad_nonCent, y = prDload)
plot(x = newDF$deltaFrq2, y = prDfreq2)
partResidFig4 <- data.frame(deltaLoad = newDF$dLoad_nonCent,
partResDeltaLoad = prDload,
deltaFreq2 = newDF$deltaFrq2,
partResDeltaFreq2 = prDfreq2)
write.csv(partResidFig4, file = "PartialResidFig4.csv", row.names = FALSE)model diagnostics
car::vif(mm7)## deltaFrq2 dLoad_nonCent
## 1.014369 1.014369par(mfrow = c(2,2))
plot(mm7, which = 1:4) # no glaring violations
par(mfrow = c(1,1))model visualization
# calculate partial residuals for deltaLoad
# these are the residuals, minus the effect of detlaLoad
y <- residuals(mm7, type = 'partial')[, "dLoad_nonCent"]
# plot partial residuals with base R plotting
plot(x = newDF$dLoad_nonCent, y = y)
partialResidDeltaLoad <- cbind(partialResidDeltaLoad, y)
colnames(partialResidDeltaLoad)[4] <- "partResDeltaMetR"
write.csv(partialResidDeltaLoad, file = "PartialResidDeltaLoad.csv", row.names = FALSE)
# visualize partial residuals (the leftmost column is what we're plotting)
scatterplotMatrix(partialResidDeltaLoad)
# double check to make sure the slope for partial residual plots are the
# same as in the original regression
summary(lm(y ~ newDF$dLoad_nonCent)) ##
## Call:
## lm(formula = y ~ newDF$dLoad_nonCent)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.17958 -0.59123 -0.06454 0.34379 3.00903
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.6463 0.7371 -2.233 0.0337 *
## newDF$dLoad_nonCent 21.4030 9.3554 2.288 0.0299 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8744 on 28 degrees of freedom
## Multiple R-squared: 0.1575, Adjusted R-squared: 0.1274
## F-statistic: 5.234 on 1 and 28 DF, p-value: 0.02991summary(mm7)##
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + dLoad_nonCent, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.17958 -0.59123 -0.06454 0.34379 3.00903
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.777e-01 7.507e-01 0.237 0.8146
## deltaFrq2 3.451e-04 6.498e-05 5.310 1.32e-05 ***
## dLoad_nonCent 2.140e+01 9.595e+00 2.231 0.0342 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8904 on 27 degrees of freedom
## Multiple R-squared: 0.5749, Adjusted R-squared: 0.5434
## F-statistic: 18.26 on 2 and 27 DF, p-value: 9.654e-06# this is what the raw data look like
plot(x = newDF$dLoad_nonCent, y = newDF$deltaMetR)
# this package plots the partial residuals
crPlot(mm7, variable = "dLoad_nonCent")
crPlot(mm7, variable = "deltaFrq2")
summary(mm7) # final model for paper##
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + dLoad_nonCent, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.17958 -0.59123 -0.06454 0.34379 3.00903
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.777e-01 7.507e-01 0.237 0.8146
## deltaFrq2 3.451e-04 6.498e-05 5.310 1.32e-05 ***
## dLoad_nonCent 2.140e+01 9.595e+00 2.231 0.0342 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8904 on 27 degrees of freedom
## Multiple R-squared: 0.5749, Adjusted R-squared: 0.5434
## F-statistic: 18.26 on 2 and 27 DF, p-value: 9.654e-06Takeaways for model 3:
- Deltafreq^2 and deltaLoad are both associated with deltaMetR.
- Holding deltafreq constant (not necessarily at 0), an increase in deltaload is associated with an increase in deltametabolicRate
- Holding deltaLoad constant (not holding load constant, and again, not holding deltaload necessarily at 0), an increase in deltaFreq^2 is associated with an increase in deltaMetabolic rate.
- We found no evidence to suggest that a change in deltaArcL^2 is associated with deltametabolicRate. (this is not saying that arcLength doesn’t affect metabolic rate).
Refref: Do regression for Average percent loading vs. delta metabolic rate / 1% load etc. It’s the last page of the document Susie sent (3 different response variables).
Covariates: avg % load, order, bee size.
(DeltaMetRate / deltaPerLoad) ~ avgPercLoad + order_1 + size_pc1
plot(deltaFrq2 ~ deltaPercLoad, data = newDF)
plot(deltaMetR ~ deltaPercLoad, data = newDF)
plot(deltaArcL2 ~ deltaPercLoad, data = newDF)
newDF <- within(newDF, {dmr_dpl = deltaMetR / deltaPercLoad})
plot(deltaMetR ~ dmr_dpl, data = newDF)
mod1 <- lm(dmr_dpl ~ avgPercLoad + order_1 + size_pc1, data = newDF)
summary(mod1)##
## Call:
## lm(formula = dmr_dpl ~ avgPercLoad + order_1 + size_pc1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.025580 -0.011566 -0.004357 0.009067 0.043391
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.1078453 0.0158826 6.790 3.31e-07 ***
## avgPercLoad -0.0010569 0.0002784 -3.796 0.000794 ***
## order_1loadedSecond -0.0131290 0.0066538 -1.973 0.059196 .
## size_pc1 0.0014674 0.0017912 0.819 0.420121
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01773 on 26 degrees of freedom
## Multiple R-squared: 0.4844, Adjusted R-squared: 0.4249
## F-statistic: 8.141 on 3 and 26 DF, p-value: 0.0005502mod2 <- update(mod1, .~. - size_pc1)
anova(mod1, mod2) # remove size_pc1## Analysis of Variance Table
##
## Model 1: dmr_dpl ~ avgPercLoad + order_1 + size_pc1
## Model 2: dmr_dpl ~ avgPercLoad + order_1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 26 0.0081723
## 2 27 0.0083833 -1 -0.00021093 0.6711 0.4201mod3 <- update(mod2, .~. - order_1)
anova(mod2, mod3) # remove order## Analysis of Variance Table
##
## Model 1: dmr_dpl ~ avgPercLoad + order_1
## Model 2: dmr_dpl ~ avgPercLoad
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 27 0.0083833
## 2 28 0.0094654 -1 -0.0010822 3.4854 0.07281 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mod3) # final model for paper##
## Call:
## lm(formula = dmr_dpl ~ avgPercLoad, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.026590 -0.011771 -0.000679 0.011114 0.046931
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.1082986 0.0158566 6.830 2.02e-07 ***
## avgPercLoad -0.0011884 0.0002735 -4.345 0.000165 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01839 on 28 degrees of freedom
## Multiple R-squared: 0.4028, Adjusted R-squared: 0.3814
## F-statistic: 18.88 on 1 and 28 DF, p-value: 0.0001654# get p-values for avgPercLoad and order
anova(mod2, update(mod3, .~. - avgPercLoad)) # p-value for avg perc load## Analysis of Variance Table
##
## Model 1: dmr_dpl ~ avgPercLoad + order_1
## Model 2: dmr_dpl ~ 1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 27 0.0083833
## 2 29 0.0158487 -2 -0.0074654 12.022 0.0001846 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# plot(I(deltaMetR) ~ avgPercLoad, data = newDF)
# plot(I(deltaMetR / avgPercLoad) ~ deltaMetR, data = newDF)
car::scatterplotMatrix(newDF[, c("avgPercLoad", "deltaMetR", "size_pc1")])
car::vif(mod1) # looks fine## avgPercLoad order_1 size_pc1
## 1.114454 1.051699 1.116654par(mfrow = c(2,2))
plot(mod3) # possible nonlinear trend in residuals
plot(mod3, which = 4)
# update model to add a non-linear term
newDF <- within(newDF, {avgPercLoad_cent = as.numeric(scale(newDF$avgPercLoad, center = TRUE, scale = FALSE))})
newDF$apl_cent2 <- with(newDF, avgPercLoad_cent^2)
m11 <- lm(dmr_dpl ~ avgPercLoad_cent + apl_cent2 + order_1 + size_pc1, data = newDF)
summary(m11)##
## Call:
## lm(formula = dmr_dpl ~ avgPercLoad_cent + apl_cent2 + order_1 +
## size_pc1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.026372 -0.011430 -0.004154 0.009170 0.044752
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.722e-02 5.369e-03 8.795 4e-09 ***
## avgPercLoad_cent -1.070e-03 2.860e-04 -3.741 0.00096 ***
## apl_cent2 7.051e-06 2.112e-05 0.334 0.74123
## order_1loadedSecond -1.373e-02 7.006e-03 -1.960 0.06125 .
## size_pc1 1.494e-03 1.824e-03 0.819 0.42065
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01804 on 25 degrees of freedom
## Multiple R-squared: 0.4866, Adjusted R-squared: 0.4045
## F-statistic: 5.925 on 4 and 25 DF, p-value: 0.0017car::vif(m11) # vif is much better with centered variables## avgPercLoad_cent apl_cent2 order_1 size_pc1
## 1.135857 1.108436 1.126036 1.118753m11a <- update(m11, .~. - size_pc1)
anova(m11, m11a) # remove size_pc1## Analysis of Variance Table
##
## Model 1: dmr_dpl ~ avgPercLoad_cent + apl_cent2 + order_1 + size_pc1
## Model 2: dmr_dpl ~ avgPercLoad_cent + apl_cent2 + order_1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 25 0.0081360
## 2 26 0.0083542 -1 -0.00021818 0.6704 0.4206m11b <- update(m11a, .~. - apl_cent2)
anova(m11a, m11b) # drop squared term## Analysis of Variance Table
##
## Model 1: dmr_dpl ~ avgPercLoad_cent + apl_cent2 + order_1
## Model 2: dmr_dpl ~ avgPercLoad_cent + order_1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 26 0.0083542
## 2 27 0.0083833 -1 -2.9042e-05 0.0904 0.7661summary(m11b)##
## Call:
## lm(formula = dmr_dpl ~ avgPercLoad_cent + order_1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.030196 -0.012515 -0.003418 0.009765 0.040799
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0474326 0.0047309 10.026 1.34e-10 ***
## avgPercLoad_cent -0.0011242 0.0002644 -4.253 0.000226 ***
## order_1loadedSecond -0.0121420 0.0065038 -1.867 0.072809 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01762 on 27 degrees of freedom
## Multiple R-squared: 0.471, Adjusted R-squared: 0.4319
## F-statistic: 12.02 on 2 and 27 DF, p-value: 0.0001846m11c <- update(m11b, .~. -order_1 )
anova(m11b, m11c) # remove order_1## Analysis of Variance Table
##
## Model 1: dmr_dpl ~ avgPercLoad_cent + order_1
## Model 2: dmr_dpl ~ avgPercLoad_cent
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 27 0.0083833
## 2 28 0.0094654 -1 -0.0010822 3.4854 0.07281 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m11c)##
## Call:
## lm(formula = dmr_dpl ~ avgPercLoad_cent, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.026590 -0.011771 -0.000679 0.011114 0.046931
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0409569 0.0033568 12.201 1.01e-12 ***
## avgPercLoad_cent -0.0011884 0.0002735 -4.345 0.000165 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01839 on 28 degrees of freedom
## Multiple R-squared: 0.4028, Adjusted R-squared: 0.3814
## F-statistic: 18.88 on 1 and 28 DF, p-value: 0.0001654## visualize model for deltametRate/avgPercLoading
par(mfrow = c(1,1))
plot(dmr_dpl ~ avgPercLoad, col = factor(order_1), data = newDF, pch = 20)
ndf <- data.frame(avgPercLoad_cent = seq(min(newDF$avgPercLoad_cent), max(newDF$avgPercLoad_cent), length.out = 200),
order_1 = as.factor(rep(c("loadedFirst", "loadedSecond"), 100)))
ndf$apl_cent2 <- ndf$avgPercLoad_cent^2
ndf$avgPercLoad <- ndf$avgPercLoad_cent + mean(newDF$avgPercLoad)
predDF <- data.frame(preds = predict(mod3, newdata = ndf), ndf)
ggplot(newDF, aes(x = avgPercLoad, y = dmr_dpl)) +
geom_point(aes(color = order_1), shape = 17) +
geom_line(data = predDF, aes(x = avgPercLoad, y = preds)) +
labs(x = "Average Load (% bodymass)",
y = "Change in Metabolic Rate (mL CO2 / hr) / \n Change in Load (% bodymass)") +
scale_color_viridis( name = "Order", discrete = TRUE, end = 0.8) +
theme(legend.position = c(0.8,0.8))
ggsave("dmr_dpl.pdf", width = 5, height = 4)(DeltaFrq^2/deltaPercLoad) ~ avgPercLoad + order_1 + size_pc1
refref: why not just put deltaPercLoad into the predictors?
newDF <- within(newDF, {df2_dpl = deltaFrq2 / deltaPercLoad})
hist(newDF$df2_dpl)
par(mfrow = c(1,1))
plot(df2_dpl ~ avgPercLoad, newDF)
plot(deltaFrq2 ~ df2_dpl, data = newDF)
mod1_f <- lm(df2_dpl ~ avgPercLoad + order_1 + size_pc1, data = newDF)
# refref: this might be better
mod1_f_DL <- lm(deltaFrq2 ~ deltaPercLoad + avgPercLoad + order_1 + size_pc1, data = newDF)
summary(mod1_f_DL)##
## Call:
## lm(formula = deltaFrq2 ~ deltaPercLoad + avgPercLoad + order_1 +
## size_pc1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2394.32 -1034.21 -43.07 1274.83 2706.53
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10068.898 3115.682 3.232 0.003438 **
## deltaPercLoad -5.827 62.476 -0.093 0.926435
## avgPercLoad -128.299 31.789 -4.036 0.000452 ***
## order_1loadedSecond -2360.025 1035.176 -2.280 0.031414 *
## size_pc1 -94.564 163.446 -0.579 0.568059
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1586 on 25 degrees of freedom
## Multiple R-squared: 0.67, Adjusted R-squared: 0.6172
## F-statistic: 12.69 on 4 and 25 DF, p-value: 8.982e-06car::vif(mod1_f_DL)## deltaPercLoad avgPercLoad order_1 size_pc1
## 3.744421 1.816786 3.182718 1.162494summary(mod1_f)##
## Call:
## lm(formula = df2_dpl ~ avgPercLoad + order_1 + size_pc1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -46.651 -21.833 -1.027 24.243 50.541
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 192.2395 27.4029 7.015 1.89e-07 ***
## avgPercLoad -2.6608 0.4803 -5.539 8.14e-06 ***
## order_1loadedSecond -36.8089 11.4800 -3.206 0.00355 **
## size_pc1 -1.5020 3.0904 -0.486 0.63103
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 30.59 on 26 degrees of freedom
## Multiple R-squared: 0.6538, Adjusted R-squared: 0.6139
## F-statistic: 16.37 on 3 and 26 DF, p-value: 3.534e-06mod2_f <- update(mod1_f, .~. - size_pc1)
anova(mod1_f, mod2_f) # remove size_pc1## Analysis of Variance Table
##
## Model 1: df2_dpl ~ avgPercLoad + order_1 + size_pc1
## Model 2: df2_dpl ~ avgPercLoad + order_1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 26 24327
## 2 27 24548 -1 -221.01 0.2362 0.631summary(mod2_f) # not final model for paper##
## Call:
## lm(formula = df2_dpl ~ avgPercLoad + order_1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -45.421 -24.372 0.984 25.355 51.493
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 188.8709 26.1342 7.227 8.99e-08 ***
## avgPercLoad -2.5918 0.4524 -5.730 4.33e-06 ***
## order_1loadedSecond -37.8192 11.1294 -3.398 0.00212 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 30.15 on 27 degrees of freedom
## Multiple R-squared: 0.6507, Adjusted R-squared: 0.6248
## F-statistic: 25.15 on 2 and 27 DF, p-value: 6.817e-07# get p-values for avgPercLoad
anova(mod2_f, update(mod2_f, .~. - order_1)) # p-value for order## Analysis of Variance Table
##
## Model 1: df2_dpl ~ avgPercLoad + order_1
## Model 2: df2_dpl ~ avgPercLoad
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 27 24548
## 2 28 35047 -1 -10499 11.547 0.00212 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1par(mfrow = c(2,2))
plot(mod2_f) # non-linearity in residuals
plot(mod2_f, which = 4)
# update model to add a non-linear term
m22 <- lm(df2_dpl ~ avgPercLoad_cent + apl_cent2 + order_1 + size_pc1, data = newDF)
summary(m22)##
## Call:
## lm(formula = df2_dpl ~ avgPercLoad_cent + apl_cent2 + order_1 +
## size_pc1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -40.78 -19.93 -4.39 11.41 50.81
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.26529 7.86521 3.975 0.000528 ***
## avgPercLoad_cent -2.84108 0.41899 -6.781 4.17e-07 ***
## apl_cent2 0.09699 0.03094 3.135 0.004358 **
## order_1loadedSecond -45.07571 10.26363 -4.392 0.000180 ***
## size_pc1 -1.13910 2.67274 -0.426 0.673614
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 26.43 on 25 degrees of freedom
## Multiple R-squared: 0.7515, Adjusted R-squared: 0.7117
## F-statistic: 18.9 on 4 and 25 DF, p-value: 2.873e-07car::vif(m22) # vif is much better with centered variables## avgPercLoad_cent apl_cent2 order_1 size_pc1
## 1.135857 1.108436 1.126036 1.118753m22a <- update(m22, .~. - size_pc1)
anova(m22, m22a) # remove size_pc1## Analysis of Variance Table
##
## Model 1: df2_dpl ~ avgPercLoad_cent + apl_cent2 + order_1 + size_pc1
## Model 2: df2_dpl ~ avgPercLoad_cent + apl_cent2 + order_1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 25 17463
## 2 26 17590 -1 -126.88 0.1816 0.6736summary(m22a)##
## Call:
## lm(formula = df2_dpl ~ avgPercLoad_cent + apl_cent2 + order_1,
## data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -42.26 -19.10 -2.71 11.28 53.26
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.61308 7.69867 4.106 0.000354 ***
## avgPercLoad_cent -2.78994 0.39507 -7.062 1.69e-07 ***
## apl_cent2 0.09756 0.03042 3.207 0.003540 **
## order_1loadedSecond -45.88915 9.92463 -4.624 9.07e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 26.01 on 26 degrees of freedom
## Multiple R-squared: 0.7497, Adjusted R-squared: 0.7208
## F-statistic: 25.96 on 3 and 26 DF, p-value: 5.55e-08par(mfrow = c(2,2))
plot(m22a) # residuals look better, though slight fan shape
plot(m22a, which = 4) ## row 22 looks highly influential
nd22 <- newDF[-22, ]
m22s <- lm(df2_dpl ~ avgPercLoad_cent + apl_cent2 + order_1, data = nd22)
summary(m22s) # no major change when we remove obs num 22##
## Call:
## lm(formula = df2_dpl ~ avgPercLoad_cent + apl_cent2 + order_1,
## data = nd22)
##
## Residuals:
## Min 1Q Median 3Q Max
## -34.038 -17.721 -6.681 11.348 54.728
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 29.57857 7.34149 4.029 0.000460 ***
## avgPercLoad_cent -2.55219 0.39112 -6.525 7.79e-07 ***
## apl_cent2 0.12123 0.03101 3.910 0.000624 ***
## order_1loadedSecond -45.07241 9.38443 -4.803 6.21e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 24.57 on 25 degrees of freedom
## Multiple R-squared: 0.7554, Adjusted R-squared: 0.7261
## F-statistic: 25.74 on 3 and 25 DF, p-value: 8.179e-08plot(m22s, which = 4)
summary(m22a) # final model for paper (though we could also report the non-centered values)##
## Call:
## lm(formula = df2_dpl ~ avgPercLoad_cent + apl_cent2 + order_1,
## data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -42.26 -19.10 -2.71 11.28 53.26
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.61308 7.69867 4.106 0.000354 ***
## avgPercLoad_cent -2.78994 0.39507 -7.062 1.69e-07 ***
## apl_cent2 0.09756 0.03042 3.207 0.003540 **
## order_1loadedSecond -45.88915 9.92463 -4.624 9.07e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 26.01 on 26 degrees of freedom
## Multiple R-squared: 0.7497, Adjusted R-squared: 0.7208
## F-statistic: 25.96 on 3 and 26 DF, p-value: 5.55e-08## visualize model for deltametRate/avgPercLoading
ndf <- data.frame(avgPercLoad_cent = seq(min(newDF$avgPercLoad_cent), max(newDF$avgPercLoad_cent), length.out = 200),
order_1 = as.factor(rep(c("loadedFirst", "loadedSecond"), 100)))
ndf$apl_cent2 <- ndf$avgPercLoad_cent^2
ndf$apl <- ndf$avgPercLoad_cent + mean(newDF$avgPercLoad)
predDF <- data.frame(preds = predict(m22a, newdata = ndf, se = TRUE), ndf)
par(mfrow = c(1,1))
plot(df2_dpl ~ avgPercLoad, col = factor(order_1), data = newDF, pch = 20)
newDF[22, ]## BeeID Mstarved S Itspan L order_H order_L M2_H
## 22 E42 0.1371 6.27e-05 0.00486 0.01097648 2 1 0.2945
## M2_L MF_H MF_L MetR_H MetR_L freq_H freq_L amp_H
## 22 0.2145 0.2761 0.2002 12.2844 11.52512 175.1276 186.0056 138.0042
## amp_L load_H load_L MT_H MT_L perLoad_H perLoad_L frce_H
## 22 117.3242 0.1482 0.07025 0.2853 0.20735 108.0963 51.23997 0.003024302
## frce_L arcL2_H arcL2_L freq2_H freq2_L U_H
## 22 0.002465804 0.0003931771 0.0002841704 30669.69 34598.09 6.945104
## U_L arcL_H arcL_L deltaLoad2 dLoad_nonCent deltaLoad
## 22 6.271125 0.01982869 0.01685735 1.0609e-06 0.07795 0.00103
## deltaFreq2Perc deltaArcL2 deltaArcL deltaFrq2 deltaMetR avgPercLoad
## 22 -69.09351 0.0001090067 0.00297134 -3928.402 0.7592781 79.66813
## deltaPercLoad size_pc1 size_pc1_2 order_1 dmr_dpl
## 22 56.85631 0.4108091 0.1687641 loadedSecond 0.01335433
## avgPercLoad_cent apl_cent2 df2_dpl
## 22 23.00355 529.1635 -69.09351ggplot(newDF, aes(x = avgPercLoad, y = df2_dpl)) +
geom_point(aes(color = order_1),shape = 17) +
geom_point(aes(size = BeeID == "E42")) + # show the influential point
geom_line(data = predDF, aes(x = apl, y = preds.fit, color = order_1)) +
# geom_line(data = predDF, aes(x = apl, y = preds.fit + 1.96*preds.se.fit, color = order_1), lty = 2) +
# geom_line(data = predDF, aes(x = apl, y = preds.fit - 1.96*preds.se.fit, color = order_1), lty = 2) +
labs(x = "Average Load (% bodymass)",
y = "Change in wingbeat freq^2 (hz^2) / \n Change in Load (% bodymass)") +
scale_color_viridis( name = "Order", discrete = TRUE, end = 0.8) +
theme(legend.position = c(0.8,0.8)) ## Warning: Using size for a discrete variable is not advised.
ggplot(newDF, aes(x = avgPercLoad, y = df2_dpl)) +
geom_point(aes(color = order_1),shape = 17) +
geom_line(data = predDF, aes(x = apl, y = preds.fit, color = order_1)) +
# geom_line(data = predDF, aes(x = apl, y = preds.fit + 1.96*preds.se.fit, color = order_1), lty = 2) +
# geom_line(data = predDF, aes(x = apl, y = preds.fit - 1.96*preds.se.fit, color = order_1), lty = 2) +
labs(x = "Average Load (% bodymass)",
y = "Change in wingbeat freq^2 (hz^2) / \n Change in Load (% bodymass)") +
scale_color_viridis( name = "Order", discrete = TRUE, end = 0.8) +
theme(legend.position = c(0.8,0.8)) 
ggsave("df2_dpl.pdf", width = 5, height = 4)(DeltaArcL^2/deltaPerLoad) ~ avgPercLoad + order_1 + size_pc1
newDF <- within(newDF, {da2_dpl = deltaArcL2 / deltaPercLoad})
plot(deltaArcL2 ~ da2_dpl, data = newDF)
par(mfrow = c(1,1))
plot(da2_dpl ~ avgPercLoad, ylab = c("deltaArcL^2 / avgPercLoad"), newDF)
mod1_a <- lm(da2_dpl ~ avgPercLoad + order_1 + size_pc1, data = newDF)
summary(mod1_a)##
## Call:
## lm(formula = da2_dpl ~ avgPercLoad + order_1 + size_pc1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.132e-06 -4.984e-07 6.244e-08 4.953e-07 1.054e-06
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.719e-06 5.312e-07 3.237 0.00329 **
## avgPercLoad 3.674e-09 9.312e-09 0.395 0.69641
## order_1loadedSecond -2.368e-07 2.226e-07 -1.064 0.29706
## size_pc1 4.140e-08 5.991e-08 0.691 0.49567
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.93e-07 on 26 degrees of freedom
## Multiple R-squared: 0.05092, Adjusted R-squared: -0.05859
## F-statistic: 0.465 on 3 and 26 DF, p-value: 0.7092mod2_a <- update(mod1_a, .~. - avgPercLoad) # I also checked for a squared term in this model (code not shown)
anova(mod1_a, mod2_a) # remove avgPercLoad## Analysis of Variance Table
##
## Model 1: da2_dpl ~ avgPercLoad + order_1 + size_pc1
## Model 2: da2_dpl ~ order_1 + size_pc1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 26 9.1429e-12
## 2 27 9.1976e-12 -1 -5.4736e-14 0.1557 0.6964summary(mod2_a)##
## Call:
## lm(formula = da2_dpl ~ order_1 + size_pc1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.103e-06 -5.079e-07 7.199e-08 4.954e-07 9.781e-07
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.919e-06 1.568e-07 12.242 1.57e-12 ***
## order_1loadedSecond -2.214e-07 2.156e-07 -1.027 0.314
## size_pc1 3.442e-08 5.634e-08 0.611 0.546
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.837e-07 on 27 degrees of freedom
## Multiple R-squared: 0.04524, Adjusted R-squared: -0.02549
## F-statistic: 0.6396 on 2 and 27 DF, p-value: 0.5353mod2_b <- update(mod2_a, .~. - size_pc1)
anova(mod2_a, mod2_b) # remove size## Analysis of Variance Table
##
## Model 1: da2_dpl ~ order_1 + size_pc1
## Model 2: da2_dpl ~ order_1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 27 9.1976e-12
## 2 28 9.3248e-12 -1 -1.2716e-13 0.3733 0.5463mod2c <- update(mod2_b, .~. - order_1)
anova(mod2c, mod2_b) # remove order## Analysis of Variance Table
##
## Model 1: da2_dpl ~ 1
## Model 2: da2_dpl ~ order_1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 29 9.6334e-12
## 2 28 9.3248e-12 1 3.0861e-13 0.9267 0.344summary(mod2c) # final mod for paper##
## Call:
## lm(formula = da2_dpl ~ 1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.206e-06 -4.613e-07 -1.433e-08 4.925e-07 9.636e-07
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.801e-06 1.052e-07 17.12 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.764e-07 on 29 degrees of freedompar(mfrow = c(2,2))
plot(mod2_a)
plot(mod2_a, which = 4)
# visualize model
ndf <- data.frame(avgPercLoad_cent = seq(min(newDF$avgPercLoad_cent), max(newDF$avgPercLoad_cent), length.out = 200),
order_1 = as.factor(rep(c("loadedFirst", "loadedSecond"), 100)))
ndf$apl_cent2 <- ndf$avgPercLoad_cent^2
ndf$avgPercLoad <- ndf$avgPercLoad_cent + mean(newDF$avgPercLoad)
predDF <- data.frame(preds = predict(mod2c, newdata = ndf, se = TRUE), ndf)
par(mfrow = c(1,1))
plot(da2_dpl ~ avgPercLoad, col = factor(order_1), data = newDF, pch = 20)
ggplot(newDF, aes(x = avgPercLoad, y = da2_dpl)) +
geom_point(aes(color = order_1), shape = 17) +
geom_line(data = predDF, aes(x = avgPercLoad, y = preds.fit)) +
# geom_line(data = predDF, aes(x = apl, y = preds.fit + 1.96*preds.se.fit, color = order_1), lty = 2) +
# geom_line(data = predDF, aes(x = apl, y = preds.fit - 1.96*preds.se.fit, color = order_1), lty = 2) +
labs(x = "Average Load (% bodymass)",
y = "Change in arc length^2 (radians^2) / \n Change in Load (% bodymass)") +
scale_color_viridis( name = "Order", discrete = TRUE, end = 0.8) +
theme(legend.position = c(0.8,0.8)) 
ggsave("da2_dpl.pdf", width = 5, height = 4)mod1_a_a <- lm(deltaArcL2 ~ deltaPercLoad + avgPercLoad + order_1 + size_pc1, data = newDF)
summary(mod1_a_a)##
## Call:
## lm(formula = deltaArcL2 ~ deltaPercLoad + avgPercLoad + order_1 +
## size_pc1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.089e-05 -2.346e-05 3.073e-06 2.598e-05 5.237e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.373e-05 6.692e-05 -0.653 0.5194
## deltaPercLoad 2.661e-06 1.342e-06 1.983 0.0584 .
## avgPercLoad -2.316e-08 6.828e-07 -0.034 0.9732
## order_1loadedSecond -2.585e-06 2.223e-05 -0.116 0.9084
## size_pc1 4.078e-06 3.510e-06 1.162 0.2564
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.405e-05 on 25 degrees of freedom
## Multiple R-squared: 0.3609, Adjusted R-squared: 0.2586
## F-statistic: 3.529 on 4 and 25 DF, p-value: 0.02046m2 <- update(mod1_a_a, .~. - avgPercLoad)
anova(mod1_a_a, m2)## Analysis of Variance Table
##
## Model 1: deltaArcL2 ~ deltaPercLoad + avgPercLoad + order_1 + size_pc1
## Model 2: deltaArcL2 ~ deltaPercLoad + order_1 + size_pc1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 25 2.8991e-08
## 2 26 2.8992e-08 -1 -1.3338e-12 0.0012 0.9732summary(m2)##
## Call:
## lm(formula = deltaArcL2 ~ deltaPercLoad + order_1 + size_pc1,
## data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.098e-05 -2.320e-05 2.827e-06 2.600e-05 5.240e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.325e-05 6.412e-05 -0.675 0.5059
## deltaPercLoad 2.633e-06 1.031e-06 2.555 0.0168 *
## order_1loadedSecond -3.028e-06 1.764e-05 -0.172 0.8650
## size_pc1 4.090e-06 3.424e-06 1.195 0.2430
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.339e-05 on 26 degrees of freedom
## Multiple R-squared: 0.3609, Adjusted R-squared: 0.2871
## F-statistic: 4.893 on 3 and 26 DF, p-value: 0.007913m3 <- update(m2, .~. - order_1)
anova(m2, m3)## Analysis of Variance Table
##
## Model 1: deltaArcL2 ~ deltaPercLoad + order_1 + size_pc1
## Model 2: deltaArcL2 ~ deltaPercLoad + size_pc1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 26 2.8992e-08
## 2 27 2.9025e-08 -1 -3.2866e-11 0.0295 0.865summary(m3)##
## Call:
## lm(formula = deltaArcL2 ~ deltaPercLoad + size_pc1, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.208e-05 -2.333e-05 1.125e-06 2.626e-05 5.136e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.183e-05 3.945e-05 -1.314 0.199952
## deltaPercLoad 2.759e-06 7.078e-07 3.898 0.000579 ***
## size_pc1 4.179e-06 3.323e-06 1.257 0.219406
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.279e-05 on 27 degrees of freedom
## Multiple R-squared: 0.3601, Adjusted R-squared: 0.3127
## F-statistic: 7.598 on 2 and 27 DF, p-value: 0.002411m4 <- update(m3, .~. - size_pc1)
anova(m3, m4)## Analysis of Variance Table
##
## Model 1: deltaArcL2 ~ deltaPercLoad + size_pc1
## Model 2: deltaArcL2 ~ deltaPercLoad
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 27 2.9025e-08
## 2 28 3.0724e-08 -1 -1.6994e-09 1.5808 0.2194summary(m4)##
## Call:
## lm(formula = deltaArcL2 ~ deltaPercLoad, data = newDF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.558e-05 -1.690e-05 -2.007e-06 2.554e-05 5.841e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.554e-05 3.764e-05 -0.944 0.35315
## deltaPercLoad 2.464e-06 6.745e-07 3.652 0.00106 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.313e-05 on 28 degrees of freedom
## Multiple R-squared: 0.3227, Adjusted R-squared: 0.2985
## F-statistic: 13.34 on 1 and 28 DF, p-value: 0.001059Session Info
sessionInfo()## R version 3.3.2 (2016-10-31)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: OS X El Capitan 10.11.6
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] plyr_1.8.4 viridis_0.3.4 visreg_2.3-0
## [4] effects_3.1-2 sjPlot_2.2.1 influence.ME_0.9-8
## [7] reshape2_1.4.2 MASS_7.3-45 gsheet_0.4.2
## [10] lme4_1.1-12 Matrix_1.2-8 car_2.1-4
## [13] ggplot2_2.2.1
##
## loaded via a namespace (and not attached):
## [1] tidyr_0.6.1 splines_3.3.2 merTools_0.3.0
## [4] modelr_0.1.0 shiny_1.0.0 assertthat_0.1
## [7] stats4_3.3.2 coin_1.1-3 yaml_2.1.14
## [10] backports_1.0.5 lattice_0.20-34 quantreg_5.29
## [13] digest_0.6.12 minqa_1.2.4 colorspace_1.3-2
## [16] sandwich_2.3-4 htmltools_0.3.5 httpuv_1.3.3
## [19] psych_1.6.12 broom_0.4.1 SparseM_1.74
## [22] haven_1.0.0 purrr_0.2.2 xtable_1.8-2
## [25] mvtnorm_1.0-5 scales_0.4.1 stringdist_0.9.4.4
## [28] MatrixModels_0.4-1 arm_1.9-3 tibble_1.2
## [31] mgcv_1.8-16 DT_0.2 TH.data_1.0-8
## [34] nnet_7.3-12 lazyeval_0.2.0 pbkrtest_0.4-6
## [37] mnormt_1.5-5 survival_2.40-1 magrittr_1.5
## [40] mime_0.5 evaluate_0.10 nlme_3.1-130
## [43] foreign_0.8-67 tools_3.3.2 multcomp_1.4-6
## [46] stringr_1.1.0 munsell_0.4.3 blme_1.0-4
## [49] grid_3.3.2 nloptr_1.0.4 htmlwidgets_0.8
## [52] labeling_0.3 rmarkdown_1.3 gtable_0.2.0
## [55] codetools_0.2-15 curl_2.3 abind_1.4-5
## [58] sjstats_0.7.1 DBI_0.5-1 sjmisc_2.2.1
## [61] R6_2.2.0 gridExtra_2.2.1 zoo_1.7-14
## [64] knitr_1.15.1 dplyr_0.5.0 rprojroot_1.2
## [67] readr_1.0.0 modeltools_0.2-21 stringi_1.1.2
## [70] parallel_3.3.2 Rcpp_0.12.9 coda_0.19-1
## [73] lmtest_0.9-34Sys.time()## [1] "2017-04-17 17:30:37 EDT"